{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 磁性质数值导数 (2)：RHF 的 GIAO 磁化率"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2020-08-30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{danger}\n",
    "\n",
    "该文档的数值导数策略有误。数值导数应当要对双电子排斥积分 (ERIs) 作改动，而非将 Fock 矩阵的贡献纳入 Core Hamiltonian 中。\n",
    "\n",
    "近日该文档将会作修订。\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算 GIAO 的 RHF 数值磁化率的程序。该文档大量参考 PySCF 的代码 [magnetizability/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/magnetizability/rhf.py) 与 [nmr/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/nmr/rhf.py)。一篇公式记号比较清晰的文章是 Laasner, Blum, et al. [^Laasner-Blum.arXiv.2018]。\n",
    "\n",
    "与上一篇文档一样，我们的讨论中所使用到的分子体系 `mol` 会是非对称的氨分子，并且取用最小基组。其 RHF 计算放在实例 `mf`，而磁化率计算实例会放在 `mf_mag`。\n",
    "\n",
    "在这篇文档中，我们仍然需要保留规范原点 (Gauge Origin) 的概念。规范原点会取在坐标原点上 `coord_orig`。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto, scf, dft\n",
    "from pyscf.prop import nmr, magnetizability\n",
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(precision=5, linewidth=150, suppress=True)\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = gto.Mole()\n",
    "mol.atom = \"\"\"\n",
    "N  0.  0.  0.\n",
    "H  0.  1.  0.2\n",
    "H  0.1 0.3 1.5\n",
    "H  0.9 0.4 -.2\n",
    "\"\"\"\n",
    "mol.basis = \"STO-3G\"\n",
    "mol.verbose = 0\n",
    "mol.build()\n",
    "coord_orig = np.zeros(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nocc, nao, nmo = mol.nelec[0], mol.nao, mol.nao"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其自洽场能量为"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-55.253540514686556"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf = scf.RHF(mol).run()\n",
    "mf.e_tot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们定义磁化率计算类为 `mf_mag`。其磁化率张量 $\\xi_{ts}$ 为："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.72497, -0.02158, -0.07514],\n",
       "       [-0.02158, -2.88877,  0.20345],\n",
       "       [-0.07514,  0.20345, -3.57485]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_mag = magnetizability.RHF(mf)\n",
    "mf_mag.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "但留意到，与上一份文档不同地，如果我们对分子作平移操作 (譬如下述平移操作是将原点移动 $(x, y, z) = (10, -10, 5) \\, \\mathsf{a.u.}$)，其磁化率仍然保持不变："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7fa032502d00>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_trans = mol.copy()\n",
    "mol_trans.set_geom_(mol.atom_coords() + np.array([[10, -10, 5]]), unit=\"AU\")\n",
    "mol_trans.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.72498, -0.0216 , -0.07515],\n",
       "       [-0.0216 , -2.8888 ,  0.20345],\n",
       "       [-0.07515,  0.20345, -3.57482]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "magnetizability.RHF(scf.RHF(mol_trans).run()).kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GIAO 规范不变原子轨道"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 基本概念"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们从上一篇文档中，得知对分子作平移操作后，磁化率可能会改变。\n",
    "\n",
    "GIAO 规范不变原子轨道 (Gauge Invariant Atomic Orbital) 是一种解决方案，它可以让分子在平移操作后，保证磁化率不变。事实上，分子旋转操作下磁化率也不变，但这会额外引入三维旋转矩阵，为了方便我们就不讨论分子旋转的情形。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GIAO 的基本思路是对于外磁场引入的微扰，其波函数更换为\n",
    "\n",
    "$$\n",
    "\\phi^\\mathrm{GIAO}_\\mu (\\boldsymbol{r}) = e^{- \\frac{i}{2} (\\boldsymbol{R}_\\mu \\times \\boldsymbol{r}) \\cdot \\boldsymbol{\\mathscr{B}}} \\phi_\\mu (\\boldsymbol{r})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，$\\phi_\\mu (\\boldsymbol{r})$ 是普通的原子轨道基组，$\\boldsymbol{R}_\\mu$ 是原子轨道 $\\mu$ 作为 Gaussian 基组的中心坐标相对于规范原点 (Gauge Origin) 向量，$\\boldsymbol{r}$ 是电子坐标，$\\boldsymbol{\\mathscr{B}}$ 是外加微扰磁场。根据其定义，我们知道，GIAO 方法必须要使用原子轨道基组下使用。需要注意，这是一个复函数轨道。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用该转换后的原子轨道就可以达到平移操作后磁化率不变的结果。关于这一点的证明，参考 Pople [^Pople-Pople.MP.1958] eq 2.5 附近的讨论。\n",
    "\n",
    "事实上，规范 (原点) **不变** 原子轨道的称呼可能是不恰当的。Pople [^Pople-Pople.DFS.1962.34] 曾在文章里使用了规范 (原点) **依赖** 原子轨道 (Gauge Dependent Atomic Orbital)；这是因为 $\\phi^\\mathrm{GIAO}_\\mu (\\boldsymbol{r})$ 实际上是依规范原点位置不同而不同的 (注意 $\\boldsymbol{R}_\\mu$ 的定义)。当然，我们仍然沿用 GIAO 的约定俗成。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 实例分析：GIAO 重叠矩阵在外磁场下的微扰"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**概念与程序准备**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们拿一个具体的例子来说明 GIAO 的计算；最简单的例子是重叠矩阵：\n",
    "\n",
    "$$\n",
    "S_{\\mu \\nu}^\\mathrm{GIAO} = \\langle \\mu | e^{- \\frac{i}{2} (\\boldsymbol{R}_{\\mu \\nu} \\times \\boldsymbol{r}) \\cdot \\boldsymbol{\\mathscr{B}}} | \\nu \\rangle\n",
    "$$\n",
    "\n",
    "其中，$\\boldsymbol{R}_{\\mu \\nu} = \\boldsymbol{R}_\\nu - \\boldsymbol{R}_\\mu$，它表示原子轨道 $\\mu$ 与 $\\nu$ 作为 Gaussian 基组的中心坐标的向量之差。之所以会出现 $- \\boldsymbol{R}_\\mu$ 的负号，是因为 $\\phi^\\mathrm{GIAO}_\\mu (\\boldsymbol{r})$ 出现在左矢时，其虚数使得在作复共轭时，指数项应当乘以负号。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们在实际计算中，只关心其一阶算符的矩阵形式：\n",
    "\n",
    "$$\n",
    "\\nabla_{\\boldsymbol{\\mathscr{B}}} S_{\\mu \\nu}^\\mathrm{GIAO} = \\langle \\mu | - \\frac{i}{2} \\boldsymbol{R}_{\\mu \\nu} \\times \\boldsymbol{r} | \\nu \\rangle\n",
    "$$\n",
    "\n",
    "其分量形式不太容易写出，但我们下面会用程序来具体地求取其分量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们先使用格点积分来计算矩阵。借助 PySCF 的 DFT 格点，我们写出下述代码：\n",
    "\n",
    "- `ni` 格点积分引擎\n",
    "\n",
    "- `grids` (50, 194) 大小的格点\n",
    "\n",
    "- `ao` $\\phi_\\mu (\\boldsymbol{r}_g)$ 即处于格点 $\\boldsymbol{r}_g$ 上 $\\phi_\\mu$ 基轨道的函数值，维度表示为 $(g, \\mu)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26896, 8)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ni = dft.numint.NumInt()\n",
    "\n",
    "grids = dft.Grids(mol)\n",
    "grids.atom_grid = (50, 194)\n",
    "grids.build()\n",
    "\n",
    "ao = ni.eval_ao(mol, grids.coords)\n",
    "ao.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**数值格点积分**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们预先已经知道，第 $\\mu = 3$ 根基轨道是第 0 个原子 (N 原子)，第 $\\nu = 6$ 跟轨道是第 2 个原子 (H 原子) (按照 0-index 计数)，那么向量 $\\boldsymbol{R}_{\\mu \\nu} = \\boldsymbol{R}_{36} = \\boldsymbol{R}_{6} - \\boldsymbol{R}_{3}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.18897, 0.56692, 2.83459])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R_36 = (mol.atom_coord(2) - coord_orig) - (mol.atom_coord(0) - coord_orig)\n",
    "R_36"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述代码中，看起来在 `coord_orig` 上有多余的代码；这是因为 $\\boldsymbol{R}_\\mu$ 本身应当是原子核坐标相对于规范原点的距离；这里我们的规范原点恰好是原点 `coord_orig`。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面我们就求取积分值 $\\nabla_{\\boldsymbol{\\mathscr{B}}} S_{36}^\\mathrm{GIAO}$："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla_{\\boldsymbol{\\mathscr{B}}} S_{36}^\\mathrm{GIAO} = \\int \\phi_3 (\\boldsymbol{r}) \\left( - \\frac{i}{2} \\boldsymbol{R}_{36} \\times \\boldsymbol{r} \\right) \\phi_6 (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r} \\simeq\n",
    "- \\frac{i}{2} \\sum_g w_g \\phi_3 (\\boldsymbol{r}_g) \\phi_6 (\\boldsymbol{r}_g) \\cdot \\boldsymbol{R}_{36} \\times \\boldsymbol{r}_g\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.+0.22657j, 0.-0.j     , 0.-0.0151j ])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- 0.5j * np.einsum(\"g, g, g, gr -> r\", grids.weights, ao[:, 3], ao[:, 6], np.cross(R_36, grids.coords))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述结果是一个纯虚数向量，其三个维度分别相当于 $(\\mathscr{B}_x, \\mathscr{B}_y, \\mathscr{B}_z)$ 的维度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，我们指出，在 PySCF/libcint 中，上述积分的虚部负值可以通过输入 `int1e_igovlp` 字符串实现。其意义可以参考 [auto_intor.cl](https://github.com/sunqm/libcint/blob/abf6948fa17e5b4ecbd26de05bf4b1d7b2b2fe3c/scripts/auto_intor.cl#L17) 与 [README](https://github.com/sunqm/libcint/blob/master/README)。留意在 libcint 中，这类积分被表示为 $i \\langle \\boldsymbol{U}_\\mathrm{g} \\mu | \\nu \\rangle$，其 $\\boldsymbol{U}_\\mathrm{g} = (U_\\mathrm{g}^x, U_\\mathrm{g}^y, U_\\mathrm{g}^z)$ 的意义相当于三维向量算符\n",
    "\n",
    "$$\n",
    "\\boldsymbol{U}_\\mathrm{g} = - \\frac{i}{2} \\boldsymbol{R}_{\\mu \\nu} \\times \\boldsymbol{r}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.22658, -0.     ,  0.01511])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol.intor(\"int1e_igovlp\")[:, 3, 6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PySCF 导出是实 (反对称) 矩阵 $i \\langle \\boldsymbol{U}_\\mathrm{g} \\mu | \\nu \\rangle$；因此在实际使用该积分时，我们就需要乘以 $- i$，让该矩阵变为复厄米矩阵，成为真正的 GIAO 变换算符一阶的矩阵形式 $\\langle \\boldsymbol{U}_\\mathrm{g} \\mu | \\nu \\rangle$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.+0.22658j, 0.+0.j     , 0.-0.01511j])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- 1j * mol.intor(\"int1e_igovlp\")[:, 3, 6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "后文我们也会用 $\\boldsymbol{U}_\\mathrm{g}$ 来简化一阶导数下的 GIAO 变换算符。\n",
    "\n",
    "当然，我们需要知道这并不是一个很严格的写法，因为 $\\boldsymbol{U}_\\mathrm{g}$ 的等式右边与 $\\mu, \\nu$ 有关；因此，使用这种简写时一定需要留意哪些轨道与该 GIAO 变换算符一阶导数有关。拿双电子积分 $( U_\\mathrm{g}^t \\mu \\lambda | \\kappa \\nu )$ 来说，由于 $U_\\mathrm{g}^t$ 作用在 $\\mu$ 上，并且 $\\lambda$ 与 $\\mu$ 在积分过程中使用相同的电子坐标，因此 $U_\\mathrm{g}^t$ 的作用对象就是 $\\mu$ 与 $\\lambda$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**从偶极积分出发给出微扰的 GIAO 重叠矩阵**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "实际上，\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\nabla_{\\boldsymbol{\\mathscr{B}}} S_{\\mu \\nu}^\\mathrm{GIAO}\n",
    "&= \\langle \\mu | - \\frac{i}{2} \\boldsymbol{R}_{\\mu \\nu} \\times \\boldsymbol{r} | \\nu \\rangle \\\\\n",
    "&= - \\frac{i}{2} \\boldsymbol{R}_{\\mu \\nu} \\times \\langle \\mu | \\boldsymbol{r} | \\nu \\rangle\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们仍然拿 $(\\mu, \\nu) = (3, 6)$ 的情形讨论。我们应当注意到，$\\langle \\mu | \\boldsymbol{r} | \\nu \\rangle$ 与偶极积分形式几乎一致 (依照不同的定义方式，偶极积分可能是该值或其相反数)，可以用字符串表示为 `int1e_r`。因此，$\\nabla_{\\boldsymbol{\\mathscr{B}}} S_{36}^\\mathrm{GIAO}$ 还可以程序化为"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.+0.22658j, 0.-0.j     , 0.-0.01511j])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- 0.5j * np.cross(R_36, mol.intor(\"int1e_r\")[:, 3, 6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**微扰的 GIAO 重叠矩阵并非“规范不变”**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们知道“规范不变”的意义是，作为最终结果的磁化率在任意的规范原点下值相同。但这并不意味着中间过程的矩阵或数值结果也相同。重叠矩阵就不满足这种“规范不变”性质："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.+0.22658j, 0.+0.j     , 0.-0.01511j])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- 1j * mol.intor(\"int1e_igovlp\")[:, 3, 6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.-0.52235j, 0.-0.65814j, 0.+0.16645j])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- 1j * mol_trans.intor(\"int1e_igovlp\")[:, 3, 6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PySCF 所使用的默认规范原点就是坐标系原点，因此我们这里都没有严格地使用 `with mol.with_common_orig(coord_orig)` 语句。当然，由于 GIAO 下的磁化率应当不受规范原点的选取变化而受到影响，因此仅从最终结果上来说，也不需要刻意地规定规范原点。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Core Hamiltonian 的程序实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 一阶 Hamiltonian Core"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "前一篇文档中，所有的一阶微扰是\n",
    "\n",
    "$$\n",
    "\\hat h {}^{(1)} (\\boldsymbol{\\mathscr{B}}) = \\frac{1}{2} \\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{r} \\times \\boldsymbol{\\hat{p}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "但当考虑到 GIAO，一阶微扰算符则应当写作\n",
    "\n",
    "$$\n",
    "\\hat h {}^{(1)} (\\boldsymbol{\\mathscr{B}}, \\boldsymbol{U}_\\mathrm{g}) | \\nu \\rangle = \\frac{1}{2} \\boldsymbol{\\mathscr{B}} \\cdot (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) \\times \\boldsymbol{\\hat{p}} + \\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{U}_\\mathrm{g} \\hat f {}^{(0)} | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "需要注意到，这里使用到了零阶 Fock 算符 $\\hat f {}^{(0)}$，该算符需要代入未受外场微扰的密度矩阵才能获得，而不是单纯地由分子构型与基组就能决定的。该项贡献的来源是 $\\langle \\Psi^{(0), \\mathrm{GIAO}} | \\hat H^{(0)} | \\Psi^{(0), \\mathrm{GIAO}} \\rangle$。\n",
    "\n",
    "因此，顺磁项对应的积分可以公式表达为\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "h_{t \\mu \\nu}^{(1)}\n",
    "&=\n",
    "\\frac{1}{2} \\langle \\mu | \\hat l_t | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}\n",
    "\\\\ & \\quad\n",
    "+ \\langle U_\\mathrm{g}^t \\mu | \\hat t | \\nu \\rangle\n",
    "+ \\langle U_\\mathrm{g}^t \\mu | \\hat v_\\mathrm{nuc} | \\nu \\rangle\n",
    "\\\\ & \\quad\n",
    "+ \\sum_{\\kappa \\lambda} ( U_\\mathrm{g}^t \\mu \\nu | \\kappa \\lambda ) D_{\\kappa \\lambda}^{(0)}\n",
    "- \\frac{1}{2} \\sum_{\\kappa \\lambda} ( U_\\mathrm{g}^t \\mu \\lambda | \\kappa \\nu ) D_{\\kappa \\lambda}^{(0)}\n",
    "- \\frac{1}{2} \\sum_{\\kappa \\lambda} ( U_\\mathrm{g}^t \\kappa \\nu | \\mu \\lambda ) D_{\\kappa \\lambda}^{(0)}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，积分字符与公式表达之间的关系为\n",
    "\n",
    "- `int1e_cg_irxp` $i \\langle \\mu | \\hat l_t | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}$\n",
    "\n",
    "- `int1e_igkin` $i \\langle \\boldsymbol{U}_\\mathrm{g} \\mu | \\hat t | \\nu \\rangle$\n",
    "\n",
    "- `int1e_ignuc` $i \\langle \\boldsymbol{U}_\\mathrm{g} \\mu | \\hat v_\\mathrm{nuc} | \\nu \\rangle$\n",
    "\n",
    "- `int2e_ig1` $i ( \\boldsymbol{U}_\\mathrm{g} \\mu \\nu | \\kappa \\lambda )$\n",
    "\n",
    "需要留意 $i$，因为这会使得上式中很多项的正负号在编写程序时是相反的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此，一阶 Core Hamiltonian `hcore_1` $h_{t \\mu \\nu}^{(1)}$ 的表达式可以用下述程序表示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm_guess = mf.make_rdm1()\n",
    "hcore_1 = 1j * (\n",
    "    - 0.5 * mol.intor(\"int1e_giao_irjxp\")\n",
    "    - mol.intor(\"int1e_igkin\")\n",
    "    - mol.intor(\"int1e_ignuc\")\n",
    "    - np.einsum(\"tuvkl, kl -> tuv\", mol.intor(\"int2e_ig1\"), dm_guess)\n",
    "    + 0.5 * np.einsum(\"tulkv, kl -> tuv\", mol.intor(\"int2e_ig1\"), dm_guess)\n",
    "    + 0.5 * np.einsum(\"tkvul, kl -> tuv\", mol.intor(\"int2e_ig1\"), dm_guess)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PySCF 中有对应的函数，生成 GIAO 下的一阶 Core Hamiltonian 矩阵 (事实上上述程序块就是从下述函数中获得的)："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(\n",
    "    hcore_1,\n",
    "    1j * nmr.rhf.make_h10(mol, dm_guess)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，我们指出，关于库伦积分，实际上应当表示为\n",
    "\n",
    "$$\n",
    "h_{t \\mu \\nu}^{(1)} \\leftarrow\n",
    "- i \\sum_{\\kappa \\lambda} ( U_\\mathrm{g}^t \\mu \\nu | \\kappa \\lambda ) D_{\\kappa \\lambda}^{(0)}\n",
    "- i \\sum_{\\kappa \\lambda} ( \\mu \\nu | U_\\mathrm{g}^t \\kappa \\lambda ) D_{\\kappa \\lambda}^{(0)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "但由于反对称性质 $( \\mu \\nu | U_\\mathrm{g}^t \\kappa \\lambda ) = - ( \\mu \\nu | U_\\mathrm{g}^t \\lambda \\kappa )$ (反映到复数的情形其实是复共轭)，导致它与对称的零阶密度矩阵 $D_{\\kappa \\lambda}^{(0)}$ 相乘并求和后必为零值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(mol.intor(\"int2e_ig1\"), - mol.intor(\"int2e_ig1\").swapaxes(-3, -4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在程序初步编写时，ERI 积分上很容易遇到角标顺序应该怎么写的问题。这是因为在实对称矩阵情形下，$\\sum_{\\kappa \\lambda} (\\mu \\kappa | \\nu \\lambda) D_{\\kappa \\lambda}$ 与 $\\sum_{\\kappa \\lambda} (\\mu \\lambda | \\kappa \\nu) D_{\\kappa \\lambda}$ 是完全相同的，但当 $\\boldsymbol{U}_g$ 作用于 ERI 后，这种性质就不满足了。一个比较容易达成正确程序编写与公式推导的技巧是，保证 $\\mu, \\kappa$ 处在 ERI 积分的复共轭位上，而 $\\nu, \\lambda$ 处在普通位上 (或者用物理的记号来讲，当需要交换 $\\langle \\mu \\kappa | \\nu \\lambda \\rangle$ 的角标时，只能在竖线左或者竖线右相互交换，不能跨线)。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二阶 Hamiltonian Core"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在前一篇文档中，所有的二阶微扰是\n",
    "\n",
    "$$\n",
    "\\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}}) = \\frac{1}{8} \\big( \\boldsymbol{\\mathscr{B}}^2 \\boldsymbol{r}^2 - (\\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{r}) (\\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{r})^\\mathrm{T} \\big)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "但当考虑到 GIAO，二阶微扰算符则应当写作\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}}, \\boldsymbol{U}_\\mathrm{g}) | \\nu \\rangle\n",
    "&=\n",
    "\\frac{1}{8} \\big( \\boldsymbol{\\mathscr{B}}^2 (\\boldsymbol{r} - \\boldsymbol{R}_\\nu)^2 - \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\big( (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) (\\boldsymbol{r} - \\boldsymbol{R}_\\nu)^\\mathrm{T} \\big) \\boldsymbol{\\mathscr{B}} | \\nu \\rangle \\\\\n",
    "& \\quad\n",
    "+ \\frac{i}{4} \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\big( (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) \\times \\boldsymbol{\\hat{p}} \\big) \\boldsymbol{U}_\\mathrm{g}^\\mathrm{T} \\boldsymbol{\\mathscr{B}} | \\nu \\rangle\n",
    "+ \\frac{i}{4}  \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\boldsymbol{U}_\\mathrm{g} \\big( (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) \\times \\boldsymbol{\\hat{p}} \\big)^\\mathrm{T} \\boldsymbol{\\mathscr{B}} | \\nu \\rangle \\\\\n",
    "& \\quad\n",
    "+ \\frac{1}{2} \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\big( \\boldsymbol{U}_\\mathrm{g} \\boldsymbol{U}_\\mathrm{g}^\\mathrm{T} \\big) \\boldsymbol{\\mathscr{B}}  \\hat f^{(0)} | \\nu \\rangle\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，$\\mathbf{T}$ 是向量算符的转置，这种转置会生成矩阵形式的磁化率。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "落实到程序中，则可以写为\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "h_{ts \\mu \\nu}^{(2)} \\cdot \\mathscr{B}_t \\mathscr{B}_s\n",
    "&= \\frac{1}{8} \\big( \\delta_{ts} \\langle \\mu | x^2 + y^2 + z^2 | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu} - \\langle \\mu | ts | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu} \\big) \\\\\n",
    "&\\quad + \\frac{1}{4} \\langle U_g^t \\mu | \\hat l_s | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu} + \\frac{1}{4} \\langle U_g^s \\mu | \\hat l_t | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu} \\\\\n",
    "&\\quad + \\frac{1}{2} \\langle U_g^t U_g^s \\mu | \\hat t | \\nu \\rangle + \\frac{1}{2} \\langle U_g^t U_g^s \\mu | \\hat v_\\mathrm{nuc} | \\nu \\rangle \\\\\n",
    "&\\quad\n",
    "+ \\frac{1}{4} \\sum_{\\kappa \\lambda} \\big( (U_g^t U_g^s \\mu \\nu | \\kappa \\lambda) + (U_g^t \\mu \\nu | U_g^s \\kappa \\lambda) \\big) D_{\\kappa \\lambda}^{(0)}\n",
    "+ \\frac{1}{4} \\sum_{\\kappa \\lambda} \\big( (U_g^t U_g^s \\kappa \\lambda | \\mu \\nu) + (U_g^t \\kappa \\lambda | U_g^s \\mu \\nu) \\big) D_{\\kappa \\lambda}^{(0)} \\\\\n",
    "&\\quad\n",
    "- \\frac{1}{8} \\sum_{\\kappa \\lambda} \\big( (U_g^t U_g^s \\mu \\lambda | \\kappa \\nu) + (U_g^t \\mu \\lambda | U_g^s \\kappa \\nu) \\big) D_{\\kappa \\lambda}^{(0)}\n",
    "- \\frac{1}{8} \\sum_{\\kappa \\lambda} \\big( (U_g^t U_g^s \\kappa \\nu | \\mu \\lambda) + (U_g^t \\kappa \\nu | U_g^s \\mu \\lambda) \\big) D_{\\kappa \\lambda}^{(0)} \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，积分字符与公式表达之间的关系为\n",
    "\n",
    "- `int1e_rr_origj` $\\langle \\mu | ts | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}$\n",
    "\n",
    "- `int1e_grjxp` $\\langle U_g^t \\mu | \\hat l_s | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}$\n",
    "\n",
    "- `int1e_ggkin` $\\langle U_g^t \\mu | \\hat t | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}$\n",
    "\n",
    "- `int1e_ggnuc` $\\langle U_g^s \\mu | \\hat v_\\mathrm{nuc} | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}$\n",
    "\n",
    "- `int2e_gg1` $(U_g^t U_g^s \\mu \\nu | \\kappa \\lambda)$\n",
    "\n",
    "- `int2e_g1g2` $(U_g^t \\mu \\nu | U_g^s \\kappa \\lambda)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "int1e_rr_origj = mol.intor(\"int1e_rr_origj\").reshape(3, 3, nao, nao)\n",
    "int1e_grjxp = mol.intor('int1e_grjxp').reshape(3, 3, nao, nao)\n",
    "int1e_ggkin = mol.intor('int1e_ggkin').reshape(3, 3, nao, nao)\n",
    "int1e_ggnuc = mol.intor('int1e_ggnuc').reshape(3, 3, nao, nao)\n",
    "int2e_gg1   = mol.intor(\"int2e_gg1\")  .reshape(3, 3, nao, nao, nao, nao)\n",
    "int2e_g1g2  = mol.intor(\"int2e_g1g2\") .reshape(3, 3, nao, nao, nao, nao)\n",
    "int2e_gg    = int2e_gg1 + int2e_g1g2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此，二阶 Core Hamiltonian `hcore_2` $h_{ts \\mu \\nu}^{(2)}$ 的表达式可以用下述程序表示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcore_2 = (\n",
    "    + 1/8 * (np.einsum(\"ts, uv -> tsuv\", np.eye(3), int1e_rr_origj.diagonal(0, 0, 1).sum(-1)) - int1e_rr_origj)\n",
    "    + 1/4 * mol.intor('int1e_grjxp').reshape(3, 3, nao, nao)\n",
    "    + 1/4 * mol.intor('int1e_grjxp').reshape(3, 3, nao, nao).swapaxes(0, 1)\n",
    "    + 1/2 * mol.intor('int1e_ggkin').reshape(3, 3, nao, nao)\n",
    "    + 1/2 * mol.intor('int1e_ggnuc').reshape(3, 3, nao, nao)\n",
    "    + 1/4 * np.einsum(\"tsuvkl, kl -> tsuv\", int2e_gg, dm_guess)\n",
    "    + 1/4 * np.einsum(\"tskluv, kl -> tsuv\", int2e_gg, dm_guess)\n",
    "    - 1/8 * np.einsum(\"tsulkv, kl -> tsuv\", int2e_gg, dm_guess)\n",
    "    - 1/8 * np.einsum(\"tskvul, kl -> tsuv\", int2e_gg, dm_guess)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 重叠矩阵的程序实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 GIAO 下，除了 Hamiltonian Core 会发生改变，重叠矩阵也同样会产生变化。\n",
    "\n",
    "一阶重叠矩阵程序中表示为\n",
    "\n",
    "$$\n",
    "S_{t \\mu \\nu}^{(1)} = \\langle U_g^t \\mu | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "二阶重叠矩阵程序中表示为\n",
    "\n",
    "$$\n",
    "S_{t \\mu \\nu}^{(2)} = \\frac{1}{2} \\langle U_g^t U_g^s \\mu | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，积分字符与公式表达之间的关系为\n",
    "\n",
    "- `int1e_igovlp` $i \\langle U_g^t \\mu | \\nu \\rangle$\n",
    "\n",
    "- `int1e_ggovlp` $\\langle U_g^t U_g^s \\mu | \\nu \\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "ovlp_1 = - 1j * mol.intor(\"int1e_igovlp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "ovlp_2 = 0.5 * mol.intor('int1e_ggovlp').reshape(3, 3, mol.nao, mol.nao)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数值导数求磁化率"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "与上一篇文档一样，我们已经获得了原子轨道形式的一阶、二阶 Core Hamiltonian、重叠积分。通过在 `scf.RHF` 的自洽场过程中重载 (override) 函数 `get_hcore` 与 `get_ovlp`，就可以得到受外磁场微扰的分子能量 $E_\\mathrm{tot} (\\mathscr{B}_x, \\mathscr{B}_y, \\mathscr{B}_z)$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm_guess = mf.make_rdm1()\n",
    "\n",
    "def hcore_mag_field(dev_xyz):\n",
    "    mf = scf.RHF(mol)\n",
    "    def get_hcore(mol_=mol):\n",
    "        hcore_total  = np.asarray(scf.rhf.get_hcore(mol_), dtype=np.complex128)\n",
    "        hcore_total += np.einsum(\"tuv, t -> uv\", hcore_1, dev_xyz)\n",
    "        hcore_total += np.einsum(\"tsuv, t, s -> uv\", hcore_2, dev_xyz, dev_xyz)\n",
    "        return hcore_total\n",
    "    def get_ovlp(mol_):\n",
    "        ovlp_total  = np.asarray(scf.rhf.get_ovlp(mol_), dtype=np.complex128)\n",
    "        ovlp_total += np.einsum(\"tuv, t -> uv\", ovlp_1, dev_xyz)\n",
    "        ovlp_total += np.einsum(\"tsuv, t, s -> uv\", ovlp_2, dev_xyz, dev_xyz)\n",
    "        return ovlp_total\n",
    "    mf.get_hcore = get_hcore\n",
    "    mf.get_ovlp  = get_ovlp\n",
    "    return mf.kernel(dm=dm_guess)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以用与上一篇文档相同的代码实现二阶梯度：\n",
    "\n",
    "$$\n",
    "\\xi_{ts} = - \\frac{\\partial^2 E_\\mathrm{tot} (\\boldsymbol{\\mathscr{B}})}{\\partial \\mathscr{B}_t \\partial \\mathscr{B}_s}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "eng_origin = hcore_mag_field((0, 0, 0))\n",
    "interval = 1e-4\n",
    "num_polar = np.zeros((3, 3))\n",
    "for t in range(3):\n",
    "    for s in range(3):\n",
    "        if t != s:\n",
    "            dev_xyzs = np.zeros((4, 3))\n",
    "            dev_xyzs[0, t] = dev_xyzs[0, s] = dev_xyzs[1, t] = dev_xyzs[2, s] =  interval\n",
    "            dev_xyzs[3, t] = dev_xyzs[3, s] = dev_xyzs[2, t] = dev_xyzs[1, s] = -interval\n",
    "            num_polar[t, s] = (\n",
    "                + hcore_mag_field(dev_xyzs[0])\n",
    "                - hcore_mag_field(dev_xyzs[1])\n",
    "                - hcore_mag_field(dev_xyzs[2])\n",
    "                + hcore_mag_field(dev_xyzs[3])\n",
    "            ) / (4 * interval**2)\n",
    "        else:\n",
    "            dev_xyzs = np.zeros((2, 3))\n",
    "            dev_xyzs[0, t], dev_xyzs[1, t] = interval, -interval\n",
    "            num_polar[t, t] = (\n",
    "                + hcore_mag_field(dev_xyzs[0])\n",
    "                + hcore_mag_field(dev_xyzs[1])\n",
    "                - eng_origin * 2\n",
    "            ) / (interval ** 2)\n",
    "num_polar *= -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.72496, -0.02158, -0.07515],\n",
       "       [-0.02158, -2.88877,  0.20344],\n",
       "       [-0.07514,  0.20345, -3.57484]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_polar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们最后与 PySCF 所给出的结果进行核验："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.72497, -0.02158, -0.07514],\n",
       "       [-0.02158, -2.88877,  0.20345],\n",
       "       [-0.07514,  0.20345, -3.57485]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_mag.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[^Laasner-Blum.arXiv.2018]: Laasner, R.; Huhn, W.; Colell, J.; Theis, T.; Yu, V.; Warren, W.; Blum, V. Molecular NMR Shieldings, J-Couplings, and Magnetizabilities from Numeric Atom-Centered Orbital Based Density-Functional Calculations. *arXiv*: [1805.12225v1](http://arxiv.org/abs/1805.12225v1).\n",
    "\n",
    "[^Pople-Pople.MP.1958]: Pople, J. A. Molecular Orbital Theory of Aromatic Ring Currents. *Mol. Phys.* **1958**, *1* (2), 175–180. doi: [10.1080/00268975800100211](https://doi.org/10.1080/00268975800100211).\n",
    "\n",
    "[^Pople-Pople.DFS.1962.34]: Pople, J. A. Nuclear Magnetic Resonance in Diamagnetic Materials. the Theory of Chemical Shifts. *Discuss. Faraday Soc.* **1962**, 34, 7. doi: [10.1039/df9623400007](https://doi.org/10.1039/df9623400007)."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
